library(tidyverse)
library(sandwich)
library(lmtest)
library(broom)
library(texreg)
library(patchwork)

options(scipen = 100)

df_high <- read_csv("Replication/data_high_random.csv") # main results, appendix E
df_mid <- read_csv("Replication/data_mid_random.csv") # appendix H
df_all <- read_csv("Replication/data_allrespondents.csv") # appendix G

#### run regressions ================================================== ####

dvar <- c(
  # perceptions of politicians
  "politician_security", "politician_diplomacy",
  "politician_welfare", "politician_edu",
  "ess_pol", "cand_gender",
  
  # political participation
  "job_pol", "ksat_politics", "follow_news", "discuss_fam", "pol_int", "pol_hard",
  
  # policy support
  "reg_hate", "enc_metoo", "reduce_military"
)
dvar_high_only <- c(
  "politician_security", "politician_diplomacy", "politician_welfare",
  "politician_edu", "ess_pol", "ksat_politics"
)

outmodel <-
  expand_grid(
    dv = dvar,
    iv = "~ co_ed",
    gender = c("Boy", "Girl"),
    sample = c("Mid Sch Random", "High Sch Random", "All")
  ) |>
  mutate(
    data = case_when(
      sample == "Mid Sch Random" ~ list(df_mid),
      sample == "High Sch Random" ~ list(df_high),
      sample == "All" ~ list(df_all)),
    data = map2(data, gender, ~.x |> filter(gender == .y)),
    iv = case_when(
      dv == "cand_gender" ~ paste(iv, "+ cand_age + cand_free + cand_home"),
      TRUE ~ iv),
    fm = paste(dv, iv, "+ (sibling_op > 0) + public + teacher_femprop_admin"),
    fm = if_else(sample == "High Sch Random", paste(fm, "+ co_ed_m"), fm)
  ) |>
  select(-iv) |>
  filter(!(sample == "Mid Sch Random" & dv %in% dvar_high_only))

outlist <- outmodel |>
  mutate(
    out_raw = case_when(
      sample %in% c("All") ~ 
        map2(fm, data, ~ lm(as.formula(.x), data = .y)),
      sample %in% c("High Sch Random", "Mid Sch Random") ~ 
        map2(fm, data, ~ lm(as.formula(.x), data = .y, weights = pull(.y, weights)))),
    out_lmtest = map(
      out_raw, ~coeftest(.x, vcov. = vcovCL, cluster = ~ school_id, save = TRUE)),
    tidied = map(out_lmtest, ~tidy(.x)),
    glanced = map(out_lmtest, ~glance(.x))
  )

#### plots for results ================================================== ####
participate <- c(
  "job_pol", "ksat_politics", "follow_news", "discuss_fam", "pol_int", "pol_hard")
politician <- c(
  "politician_security", "politician_diplomacy", "politician_welfare", "politician_edu")
politician_agg <- c("ess_pol", "cand_gender")
policies <- c("reg_hate", "enc_metoo", "reduce_military")
dvfct <- c(participate, politician, politician_agg, policies)

outlist_plot <- outlist |>
  mutate(
    coefstd = map(tidied, ~.x |> filter(term == "co_edTRUE") |> select(estimate, std.error))
  ) |>
  unnest_wider(coefstd) |>
  mutate(
    lower = estimate - qnorm(0.95)*std.error,
    upper = estimate + qnorm(0.95)*std.error,
    dv = factor(dv, levels = dvfct),
    dv_policy = case_when(
      dv == "reg_hate" ~ "Support for Policy\nPromoting Amicable Relations",
      gender == "Boy" & dv == "enc_metoo" ~ "Support for Policy\nBenefiting Outgroup",
      gender == "Girl" & dv == "reduce_military" ~ "Support for Policy\nBenefiting Outgroup",
      gender == "Boy" & dv == "reduce_military" ~ "Support for Policy\nBenefiting Ingroup",
      gender == "Girl" & dv == "enc_metoo" ~ "Support for Policy\nBenefiting Ingroup"),
    gender = case_when(
      gender == "Boy" ~ "Boys",
      gender == "Girl" ~ "Girls"),
    gender = factor(gender, levels = c("Boys", "Girls"))
  )

varlabs <- c(
  "follow_news" = "Follow News",
  "discuss_fam" = "Discuss\nSocial Issues\nwith Family",
  "pol_int" = "Political\nInterest",
  "pol_hard" = "Political\nEfficacy",
  "job_pol" = "Consider\nBeing\na Politician",
  "ksat_politics" = "Taking\nPolitics\nin KSAT",
  
  "ess_pol" = "Gendered Perceptions\nof Competence",  
  "politician_diplomacy" = "Men Politicians\nBetter at Diplomacy",
  "politician_security" = "Men Politicians\nBetter at Security",
  "politician_welfare" = "Women Politicians\nBetter at Welfare",
  "politician_edu" = "Women Politicians\nBetter at Education",
  "cand_gender.reverse" = "Prefer \n a Man Politician",
  
  "enc_metoo" = "Encourage\n#MeToo",
  "reduce_military" = "Reduce\nMale\nMilitary Service",
  "reg_hate" = "Regulate\nGender-based\nHate Speech"
)

mytheme <- 
  theme(
    plot.margin = margin(2,2,2,2),
    strip.text = element_text(size = 12, face = "bold", margin = margin(3, 3, 3, 3)),
    axis.text.x = element_text(size = 12),
    axis.title = element_text(size = 12),
    panel.background = element_rect(fill = NA, linewidth = 0.1),
    panel.grid.major = element_line(color = "gray", linewidth = 0.1),
    legend.position = "bottom", legend.key = element_rect(fill = NA, color = NA)
  )
theme_set(mytheme)

mycolors <- c(Boys = "#0057B7FF", Girls = "#DD5129FF")
myshapes <- c(Boys = 19, Girls = 17)

plotingr <- function(){
  list(
    geom_hline(yintercept = 0, linetype = "solid",
               color = "black", linewidth = 0.5, alpha = 0.5),
    geom_pointrange(
      aes(y = estimate, ymin = lower, ymax = upper, x = dv, color = gender, shape = gender),
      size = 1, alpha = 0.9, position = position_dodge(width = 0.3)),
    scale_color_manual(values = mycolors),
    scale_shape_manual(values = myshapes),
    scale_x_discrete(labels = varlabs, name = NULL),
    labs(y = expression(paste("Effect of Co-ed (baseline: Single-gender School)", sep = "")),
         color = NULL, shape = NULL)
  )
}

## Fig 1 -----
ggplot(filter(outlist_plot, dv %in% politician, sample == "High Sch Random"))+
  plotingr()+
  ylim(-0.6, 0.7) # to fix the y axis for the comparison between results using 'High Sch Random' sample and 'All' sample

## Fig 2 -----
ggplot(filter(outlist_plot, dv %in% politician_agg, sample == "High Sch Random"))+
  plotingr()+
  ylim(-1.9, 1.9)

## Fig 3 -----
ggplot(filter(outlist_plot, dv %in% participate, sample == "High Sch Random"))+
  plotingr()+
  theme(axis.title.y = element_text(angle = 90, size = rel(0.95)))+
  ylim(-1.3, 1.3)

## Fig 4 -----
reg <- ggplot(filter(outlist_plot, dv == "reg_hate", sample == "High Sch Random"))+
  plotingr()+
  facet_wrap(vars(dv_policy), nrow = 1)+
  theme(
    axis.text.x = element_text(size = 14), axis.title = element_text(size = 14),
    strip.text = element_text(size = 12, face = "bold", margin = margin(3, 3, 3, 3)))+
  ylim(-1.6, 1.6)

outg <- ggplot(filter(outlist_plot, dv_policy == "Support for Policy\nBenefiting Outgroup",
                      sample == "High Sch Random"))+
  plotingr()+
  facet_wrap(vars(dv_policy), nrow = 1)+
  theme(
    axis.text.x = element_text(size = 14), axis.title = element_text(size = 14),
    strip.text = element_text(size = 12, face = "bold", margin = margin(3, 3, 3, 3)))+
  ylim(-1.6, 1.6)

reg+outg+plot_layout(axis_titles="collect")

#### tables for results ================================================== ####

mycoef = list(
  'I(gender == "Boy")TRUE:co_edTRUE' = "Co-ed School:Boy",
  "co_edTRUE" = "Co-ed School",
  "co_ed_mTRUE" = "Co-ed Mid School",
  'I(gender == "Boy")TRUE' = "Boy",
  "sibling_op > 0TRUE" = "Any Siblings with Different Genders = 1",
  "publicTRUE" = "Public School",
  "teacher_femprop_admin" = "Proportion of Women Teachers",
  "cand_age60대" = "Candidate Age = 60s",
  "cand_free수학교육강화" = "Candidate Policy = Improve Math Curriculum",
  "cand_home화천" = "Candidate Hometown = Hwa-cheon",
  "(Intercept)" = "Intercept"
)

mytable <- function(df, headernm = NULL, modelnm){
  n <- df |> pull(out_raw) |> map_dbl(~.x$model |> nrow())
  r2 <- df |> pull(glanced) |> map_dbl(~.x$r.squared)
  gof <- list("Num. of Obs." = n, "R Squared" = r2)
  
  texreg(
    df |> pull(out_lmtest),
    custom.gof.rows = gof,
    stars = c(0.01, 0.05, 0.1),
    custom.header = headernm,
    custom.model.names = modelnm,
    custom.coef.map = mycoef
  )
}

## appendix E.1 -----
tb_politician <- outlist_plot |>
  filter(dv %in% politician, sample == "High Sch Random") |>
  arrange(dv)

mytable(tb_politician,
        headernm = list("Men Politicians Better at Security" = 1:2,
                        "Men Politicians Better at Diplomacy" = 3:4,
                        "Women Politicians Better at Welfare" = 5:6,
                        "Women Politicians Better at Education" = 7:8),
        modelnm = rep(c("Boys", "Girls"), 4))

## appendix E.2 -----
tb_politician_agg <- outlist_plot |>
  filter(dv %in% politician_agg, sample == "High Sch Random")

mytable(tb_politician_agg,
        headernm = list("Gendered Perception of Competence" = 1:2,
                        "Prefer a Man Politician" = 3:4),
        modelnm = rep(c("Boys", "Girls"), 2))

## appendix E.3 -----
tb_participate <- outlist_plot |>
  filter(dv %in% participate, sample == "High Sch Random", gender == "Boys") |>
  arrange(dv)

mytable(
  tb_participate,
  modelnm = c("Would Consider Being a Politician",
              "Taking Politics in KSAT",
              "Follow News",
              "Discuss Social Issues with Family", 
              "Pol. Interest", "Pol. Efficacy"))

tb_participate <- outlist_plot |>
  filter(dv %in% participate, sample == "High Sch Random", gender == "Girls") |>
  arrange(dv)

mytable(
  tb_participate,
  modelnm = c("Would Consider Being a Politician",
              "Taking Politics in KSAT",
              "Follow News",
              "Discuss Social Issues with Family",
              "Pol. Interest", "Pol. Efficacy"))

## appendix E.4 -----
tb_policies <- outlist_plot |> filter(dv %in% policies, sample == "High Sch Random") |>
  arrange(dv)

mytable(
  df = tb_policies,
  headernm = list("Regulate Gender-based Hate Speech" = 1:2,
                  "Encourage #MeToo" = 3:4,
                  "Reduce Male Military Service" = 5:6),
  modelnm = rep(c("Boys", "Girls"), 3))
